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Abstract 

We present a simulation of the SU(3) spin model with chemical 
potential using a recently proposed flux representation. In this rep- 
resentation the complex phase problem is avoided and a Monte Carlo 
simulation in terms of the fluxes becomes possible. We explore the 
phase diagram of the model as a function of temperature and chemical 
potential. 
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1 Introductory remarks 



The complex phase problem (fermion sign problem) is the central obstacle 
that has held up Monte Carlo simulations of finite density lattice QCD for 
two decades. When a chemical potential is coupled the fermion determinant 
becomes complex and cannot be directly used as a probability weight in a 
Monte Carlo simulation. The alternative approaches that were explored, 
such as reweighting, power series expansion, strong coupling/large mass ex- 
pansion or analytic continuation from imaginary chemical have had only 
limited success so far. 

For several systems which are simpler than QCD the complex phase 
problem was solved by rewriting the theory in terms of new degrees of free- 
dom [TJ [2] where the complex phase problem is avoided and simulations 
with new techniques such as worm algorithms [3] become possible. Several 
of these models were also explored with another method for systems with 
sign problems, the complex Langevin approach [4J [H (H [7] . 

In this article we present results for the SU(3) spin model [3], which at 
vanishing external magnetic field is expected to describe the deconfinement 
transition [8] of pure SU(3) gauge theory. Adding the magnetic field terms 
takes into account the leading corrections from the fermion determinant and 
couples the chemical potential (see [9] for a treatment of the model beyond 
these leading contributions). The model has been studied successfully us- 
ing the complex Langevin approach [U [5j [7] and has a flux representation 
[2] where the complex phase problem is avoided. In this paper we discuss 
a Monte Carlo simulation of the SU(3) spin model based on the flux rep- 
resentation and explore the phase diagram of the model as a function of 
temperature and chemical potential. 



2 The model and its flux representation 



The action of the SU(3) spin model is given by 
S = -J2 \ T ^2[P( x ) p (x + z>)* + c.c. + k e^P(x) + e~^P(x] 

x \ u=l / 

The degrees of freedom P(x), which we sometimes will refer to as Polyakov 
loops, are the traced SU(3) variables P{x) = Tr L(x) with L(x) € SU(3). 
They are attached to the sites x of a three-dimensional cubic lattice which we 
consider to be finite with periodic boundary conditions. By i> we denote the 
unit vector in ^-direction, with u = 1, 2, 3. The first two terms are a nearest 
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neighbor interaction between adjacent Polyakov loops. The parameter r 
depends on the temperature (it increases with temperature) and is real and 
positive. The real and positive parameter k is proportional to the number of 
flavors and depends on the fermion mass (it decreases with m q ) and (i is the 
chemical potential. For later use we introduce the abbreviations rj = ne^ 
and rj = Ke _/ \ 

The grand canonical partition function of the model described by ([I]) is 
obtained by integrating the Boltzmann factor e~ s ^ over all configurations 
of the Polyakov loop variables. The corresponding measure is a product over 
the reduced Haar measures dP(x) at the sites x. Thus 



Without the magnetic term, i.e., for k = 0, the system has a low tem- 
perature phase where the expectation value (P) for the spatially summed 
Polyakov loop P = Y2 X P( x ) vanishes which is interpreted as the confined 
phase of QCD. At t c ~ 0.137 the system undergoes a first order deconfine- 
ment transition into a phase which is characterized by (P) ^ (deconfined 
phase). For small k and = the first order line persists ending in a sec- 
ond order endpoint [4j. We will show here that for fi > the first order 
transition is weakened further, i.e., the endpoint shifts towards smaller k. 

Applying high temperature expansion techniques, the partition function 
can be rewritten in terms of new degrees of freedom, the flux variables. The 
general steps to obtain the flux representation are [2]: 

• The Boltzmann factor is written as a product over all nearest neighbor 
terms and a product over all sites for the magnetic terms, respectively. 
Each individual exponential is then expanded in a power series. 

• The emerging products can be reorganized, such that at each lattice 
site x one has a (reduced) Haar measure integral over powers of P(x) 
and P(x)* . The exponents are combinations of the expansion variables 
of the individual exponential terms of the Boltzmann factors. 

• Based on techniques presented in [10] these integrals (moments of the 
one-link integrals) can be solved in closed form. The integrals give 
rise to local constraints for the allowed combinations of the expansion 
coefficients of the exponentials. 

Obviously these steps are a straightforward application of textbook high 
temperature expansion techniques applied to a somewhat unusual spin sys- 
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tern. The final result for the flux representation of the partition sum is given 
by 0: 



l x J i x v \ 




{1,1} {<>,*} 



In this form the partition sum is a sum over all configurations of the two 
sets of dimer variables l x ,ujx,u £ [0, +oo) , living on the links (x, v) and the 
monomers s x ,s x G [0, +oo), living on the sites x. The first two products in 
([5]) are weight factors that come from the expansion of the individual expo- 
nentials and it is obvious that in the second product the chemical potential 
ji enters via the powers of r\ = ne^ 1 and rj = ne~^. 

The last product is over the group integrals at each lattice site x, 

I(n,n) = [ dL (TrL) n (Trtff , (4) 
JSU{3) 

where dP denotes SU(3) Haar measure and n and n are non-negative inte- 
gers. These integrals can be evaluated in closed form [2j and turn out to be 
real and non-negative. However, the I(n,n) are non-zero only if the triality 
condition (n — n) mod 3 = is obeyed. In the flux representation ([3]) for the 
partition sum the arguments of the I are the summed fluxes f x and f x at 
the sites x of the lattice defined by 

3 3 

fx — ^ [ lx,u ~t~ ~^x— i>,v \ ~{~ Sx j fx — ^ [lx,u ~{~ ^x—v,v \ ~{~ &x ■ (5) 

v=l v=\ 

The triality condition for non-vanishing weights I then reads 

(fx - 7x)mod3 = 0, (6) 

which introduces a constraint for the allowed values of the dimer and monomer 
variables at each lattice site x. 

Thus in the form ([3|) the partition sum is a sum over configurations of 
the dimers l x ,ujx,u £ [0, +oo) , and the monomers s x ,s x € [0, +oo) with 
weight factors that are real and non-negative. For a successful solution to 
the complex phase problem one still has to find a Monte Carlo update which 
only generates configurations that obey the constraint ([6]). Otherwise the 
sample will be dominated by configurations with zero weight, and the signal 
will vanish exponentially with the volume. 
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3 Reparametrization and observables 



A first version of our Monte Carlo update of the flux representation was set 
up directly in terms of the flux representation given in ([3|) . It turned out that 
this algorithm did not perform very well, and we found that a reparametriza- 
tion of the dimer and monomer variables gives rise to an alternative flux 
representation which allows for a considerably better algorithm. 

We introduce new integer valued dimer variables k X;V E [0, +00) and 
k x ,v € (—00, +00), which are related to the old dimer variables via l X}V — 
Ixv = k x u and l xv + l xv = \k xv \ + 2k xv . The various terms in the partition 
sum change according to 

Similarly we define new integer valued monomer variables r x £ [0, +00) 
and r x € (—00, +00) which are related to the old monomer variables via 
s x — s x = r x and s x + s x = \r x \ + 2r x . The corresponding replacements in 
the partition sum are: 

e tJ,(s x -s x ) _^ e^ Tx , 

•Sx' S x l y (I '"a: I "l - ^)'^' • 

Also the fluxes f x and f x may be rewritten in terms of the new dimer and 
monomer variables 



/. = £ 

V 

7, = E 



I k Xj i/ 1 -|- 1 k x — v ^y I i j , , k xv k x ... 



kx.v "I - k x — v v + 



I k x ,u I I k x — 1 j 7 k X p kx—i/^i, 

— |- K x v + K-x—vm ^ 



\r x \+r x 
H ^ h , 



+ r ^ r rx + r x . 



(7) 



After the reparametrization the partition function thus reads: 



n 



(8) 

where the sum is now over all configurations of the new dimer variables 
k x ,u € [0, +00) and k x>v € (—00, +00) and the new monomer variables 
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fx € [0, +00) and r x £ (—00, +00). For the fluxes f x and f x now the 
expressions ([7]) are used. 

It is interesting to note that in terms of the new variables the triality 
constraint (f x — / x )mod3 = from Eq. ([6]) turns into 

( y^JFw ~ k x -v,v] + rx) mod 3 = 0. (9) 

Obviously the reparametrization simplified the constraint and only the vari- 
ables k x ,v and r x enter the constraint. The variables k x ,u and r x can be 
varied freely without generating zero weight configurations. We stress that 
also the chemical potential only couples to the r x monomers. We attribute 
the aforementioned better performance of the Monte Carlo simulation after 
the reparametrization to both these properties: The fact that the constraint 
restricts only the bared variables k x , v and r x and that the chemical potential 
couples only to the r x . 

The observables we consider are obtained as derivatives of In Z with 
respect to the physical parameters r, /i and k, and also with respect to 
the combined parameters rj = ne^ and rj = ne~^. Defining the following 
abbreviations for the sums of dimer and monomer variables, 

k = ^yk x+v \ + 2k x+u \ , r = J2^x , r= J2 f * ' 1^1 = Yl i f *i ' ( 10 ) 

X,V XX X 

one easily checks the following auxiliary identities for the first and second 
derivatives of the partition sum in the flux representation 



(fflf^*)). (n, 



Id 1 , 

z7Tt Z = -t IK 




Zdrj 2 
Z drf rf 
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From the original form (JTJ) of the action one may identify the physical in- 
terpretation of the observables obtained by derivatives of In Z. Using the 
auxiliary identities (jlip . one obtains the following flux representations for 
the internal energy U, the heat capacity C, the Polyakov loop expectation 
value (P), and the Polyakov loop susceptibility \Pi 



U = — 



c = — 



(P) 



XP 



d_ d_ _d_ 

dr ^ drj ^ drj 

r* — + 2 ^-+~ 2 ^- 
dr 2 drj 2 drf 



Z=(K+\R\ + 2R) 



-2T7J- 



o 2 



■+2rjrj- 



d 2 



drdrf drjdfj 
[(K + \R\ + 2R) - U] 2 - (K + |1| + 2R) 
R\ + R 



-2rrj 



d 2 



drdfj 



Z-U 2 



-—Z = - 

Z drj rj 



Z drj'- 



■Z-{Pf 



+ R 

U\R\+R 
rj 



+ R 



(P) 



1 f\R\+R 



+ R) 



Obviously all our observables can be computed as expectation values of sums 
of dimers and monomers and their second moments. 



4 MC algorithm and setup of the simulation 

After the reparametrization in the previous section the partition sum is a 
sum over configurations of two classes of variables: The "unbared" dimer 
and monomer variables k XtU ,r x € [0, +00) and the "bared" dimers and 
monomers, k XjU , r x £ (— 00, +00). 

The unbared variables k x>u ,r x are not subject to any constraint and we 
simply update them by randomly raising or lowering them by 1 and accept- 
ing this change with the usual Metropolis probability. Applying this step 
once to all /c Xj „ constitutes one sweep for the unbared dimers and similarly 
one sweep for the unbared monomers is to run through all r x . 

The bared variables k x „, r x must obey the constraint ([9]). The constraint 
forces the combined flux of bared dimer and monomer variables at every 
lattice point to be a multiple of 3. A starting configuration that obeys the 
constraint is given by setting all bared dimers and monomers to zero. The 
following 5 update steps leave the constraints intact and can be seen to 
give rise to an ergodic update. Each individual step is accepted with the 
Metropolis probability. 
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1. The value of a dimer variable k xv is randomly changed by ±3. For the 
sites x and x + v the right hand side of changes by three units and 
thus the constraints at the two sites remain intact. Again full sweeps 
through all k Xjl/ are implemented. 

2. The value of a monomer variable r x is randomly changed by ±3. For 
the site x the right hand side of @ changes by three units and thus the 
constraint at x remains intact. Full sweeps through all r x are done. 

3. Along an oriented plaquette all dimer variables k Xjl , are changed by 
±1 according to their orientation in the plaquette. For all corners of 
the plaquette the left hand side of ([9]) remains unchanged. A sweep is 
defined as offering this step to all plaquettes on the lattice. 

4. The value of a k xu is randomly changed by ±1 and the monomers at 
its endpoints are changed accordingly: r x r x T 1, ^x+u — > ± !■ 
The left hand sides of ([9]) remain unchanged at the two endpoints of 
the dimer. Again we define a sweep by offering this change to all bared 
dimers and the monomers at their endpoints. 

5. All values of dimer variables k x>u that sit on a straight loop which 
closes around the periodic boundary are randomly changed by ±1 
(same value for all of them). For all sites on the loop the constraint 
remains unchanged. A sweep is defined as offering this change to all 
loops in the 3 directions. At k = this step is necessary for ergodicity. 

A combined full update sweep consists of a sweep for the unbared dimers, 
a sweep for the unbared monomers, and one of each of the 5 update sweeps 
for the bared variables. 

Our simulations were done on lattices with sizes between 6 3 and 20 3 . 
We typically used statistics of up to 10 6 configurations, separated by 10 full 
update sweeps for decorrelation and 10 6 sweeps for equilibration. All errors 
we quote are statistical errors determined with the Jackknife method. 

5 Comparison to other approaches 

As a first test of our numerical results we compare the outcome of the flux 
simulation at fi = to the results from a conventional simulation in the 
spin representation (which at /i = is possible). In Fig. Q] we study on 
two different volumes the observables U,C, (P) and \Pi normalized by the 
number of sites V, as a function of the temperature r at k = 0.005, /i = 0. 
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Figure 1: Comparison of results from a simulation in the flux representation 
(filled symbols) to data from a simulation in the conventional approach 
(empty symbols). We show the results for U, C, (P) and xp as a function of 
r at k = 0.005 and fi = and compare two volumes. The observables are 
normalized by the number of lattice points V . 



For all four observables we find good agreement of the results from the 
flux simulation and the data from the conventional approach. This confirms 
that the mapping from the conventional representation to the flux degrees 
of freedom is correct, that the observables are properly represented in the 
flux picture and that the flux Monte Carlo simulation works for jj, = 0. To 
check our approach and its implementation also at \i ^ 0, we compare our 
results also to the outcome of a perturbative r-expansion. This comparison 
will be discussed in the next section. 

Another interesting test of our approach is a comparison to the recently 
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Figure 2: Comparison of (P + P*) /2V from the flux simulation (filled sym- 
bols) to the results from the complex Langevin approach (empty symbols 
and two high accuracy data points are marked with crosses). For \i = we 
also added the results from a simulation in the conventional spin approach 
(asterisks). We compare data at different values of r as a function of fi 2 for 
k = 0.02 on lattices of size 10 3 . 

published results from a complex Langevin simulation of the SU(3) spin 
model [7] (for older results with the complex Langevin approach see [HE]). 
Such a comparison can be done also for non-zero chemical potential and thus 
tests our approach as a function of all three couplings, r, k and fi. In turn, 
an agreement of the complex Langevin and the flux results is an important 
test also for the complex Langevin approach which currently sees a lot of 
interesting development. 

Fig. [2] shows that the flux results and the data from complex Langevin 





\i = 1.0 


H = 1.0 


H = 3.0 


fi = 3.0 




(P) 


(P*) 


(P) 


(P*) 


complex Langevin 
flux simulation 


0.2419(19) 
0.2416(11) 


0.3605(13) 
0.3604(13) 


1.70615(27) 
1.70627(14) 


1.74590(24) 
1.74683(17) 



Table 1: Comparison to complex Langevin results with stepsize extrapola- 
tion [7J for V = 10 3 , r = 0.125, k = 0.02 and \i as listed in the table. 
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agree very well for a wide range of parameters, and for fi = also with 
the results from a conventional simulation. A slight discrepancy is seen 
for the r = 0.132 data when the complex Langevin calculation is done 
using the lowest-order Euler discretization at a single stepsize e = 0.0005, 
without performing a zero-stepsize extrapolation. Since observables in this 
discretization depend linearly on the stepsize, the authors of [7J provided us 
with improved data points for two values, /u 2 = 0.0 and 0.2 at r = 0.132, 
where the higher-order algorithm discussed in [7J was used, in which stepsize 
corrections are essentially absent. These data points are represented by 
crosses in the plot and they nicely match the flux simulation results and for 
fi = also the data from the conventional approachQ 



6 Phase diagram 

Let us now come to the analysis of the phase diagram as a function of r and 
\x. We first present the results in Figs. [3] and [Hand discuss the details of its 
determination subsequently. In Fig. [3] we show the positions of the maxima 
of xp m the r-fj, plane (symbols connected by dotted lines) for four different 
values of k. We used lattices of size 10 3 or a combination of 10 3 , 16 3 and 20 3 
lattices for the critical points and additional checks at some of the parameter 
values. The solid curves at the bottom are the positions of the maxima from 
a perturbative expansion in r and obviously the Monte Carlo data nicely 
approach these results. This again illustrates the correctness and accuracy 
of our simulation in the flux approach. 

The horizontal dashed line indicates the critical value r c ~ 0.137 of 
the theory without external field, i.e., at k = 0. At this temperature the 
k = model undergoes a first order phase transition from a confining phase 
({P) = for r < t c ) into a deconfined phase with (P) ^ 0. For [i = the 
first order transition persists for sufficiently small k. Above some k c , which 
defines a critical endpoint, only a crossover type of behaviour is observed. 
We estimated this endpoint to be at k c = 0.016(2), r c = 0.1331(1). It is 
marked with a red cross in Fig. [3l For values of k smaller than k c a first 
order transition separates the confined and the deconfining phase at some 
critical r c € [0.133, 0.137] which depends on k. 

When the chemical potential is turned on, the transition is weakened 
further, i.e., the position of the endpoint is shifted towards smaller values 
of k. We determined the endpoint for one other line in the phase diagram, 

We thank Gert Aarts and Frank James for providing us with the data of [7] and 
extensive communication on their complex Langevin simulation of the SU(3) spin model. 
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Figure 3: The phase diagram as a function of r and /i for different values 
of k. The phase boundaries (symbols connected with dotted lines) were 
determined from the maxima of \P on 10 3 lattices, or from a combined 
analysis on 10 3 , 16 3 and 20 3 lattices. The solid curves at small r are the 
results of an expansion in r. The horizontal dashed line marks the position 
of the critical r for the k = case, and the two crosses indicate the positions 
of the endpoints for the n = and k = 0.005 cases (see the text for details). 
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Figure 4: Comparison of the maxima of xp (triangles) and C (diamonds). 
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the curve for n = 0.005. For this curve we find the endpoint to be at r c = 
0.134(1), fj, c = 1.53(10). Also this endpoint is marked with a red cross in 
Fig. [3l Points on the k = 0.005 curve to the left of that endpoint correspond 
to a first order transition, while points to the right are characterized by 
crossover behavior. 

In Fig. [J] we compare the positions of the maxima of xp (triangles) to the 
positions of the maxima of the heat capacity C (diamonds). The comparison 
is done for three values of k. For parameter values where we determined 
first order transitions (upper left corner) the positions of the maxima agree 
for the two second derivatives of the free energy. However, once we enter the 
region of crossover behavior without any singularities there is no reason to 
expect that the maxima coincide. Indeed we observe this scenario: As one 
increases /i, the curves for the maxima of xp and C start to separate and 
the right half of the phase diagram has only a very broad crossover between 
the (P) ~ region (left of the curves) and the (P) ^ phase. 

Let us now analyze the behavior of our observables near the transi- 
tion/crossover lines in more detail. In Fig. [5] we show the Polyakov loop 
susceptibility xp and the heat capacity C as a function of r for k = 0.015 
and fi = 0. This value of k is in the range where for fi = one still observes 
a first order transition for some critical value r c . The figure shows that both, 
Xp and C develop a sharp peak which scales with the volume. Furthermore, 
both these second derivatives of the free energy peak at the same position, 
and we can read off r c ~ 0.1333. 

The situation is different when we look at a larger value of k. In Fig. [6] 
we show again xp and C as a function of r, but now at re = 0.04 and fi = 0. 
Although the two quantities still show a maximum, it is obvious that the 
height of the peak does not scale with the volume, and the curves for our 
three volumes fall on top of each other. Thus we conclude that at k = 0.04 
we only observe a crossover type of transition. We furthermore point out, 
that the positions of the maxima of xp and C do not coincide, another 
indication for a crossover type of behavior. 

We cross-checked the different behavior (first order versus crossover) at 
k = 0.015 and 0.040 by inspection of histograms for the distribution of the 
internal energy U. In Fig. [7] we show the distribution of U normalized by the 
volume, for the n = 0.015, fi = case that corresponds to Fig. We use two 
values of r close to the critical value r c ~ 0.1333 that we determined from 
the position of the peaks of xp and C. The familiar two-state signal that 
characterizes first order transitions is clearly established. The situation is 
different for k = 0.040, /i = (Fig. [8]). Here we observe only a single 
peak in the histogram which changes smoothly with varying r, a behavior 
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Figure 5: Polyakov loop susceptibility and heat capacity as a function of r 
for k = 0.015 and jjl = 0. We compare three different volumes. 
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Figure 6: Polyakov loop susceptibility and heat capacity as a function of r 
for k = 0.040 and /i = 0. We compare three different volumes. 
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k = 0.015 - x = 0.1333 - (i = 0.0 K = 0.015 - T = 0.1334 - n = 0.0 




Figure 7: Histograms for the distribution of the values for the internal energy 
U /V at k = 0.015 and fi = for two different values of r close to the critical 
value. 
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Figure 8: Histograms for the distribution of the values for the internal energy 
U /V at k = 0.040 and [i = for two different values of r close to the position 
of the maximum of C. 



14 




Figure 9: Polyakov loop susceptibility and heat capacity as a function of /i 
for k = 0.02 and r = 0.12. We compare three different volumes. 



characteristic of a crossover. 

Thus from the absence of volume scaling (Fig. [5]) and the single peak 
histogram (Fig. [8} we conclude that for k = 0.04 we do not observe first 
order behavior. This finding is different from a mean field analysis of the 
SU(3) spin model [IT] where a first order behavior was reported for values 
of k up to a critical endpoint at k = 0.059. Our results show clearly that the 
endpoint must be below k = 0.04 and our best estimate is the value given 
already above in the discussion of Fig. E2 k c = 0.016(2). 

We conclude our discussion of the phase diagram with an example at 
different values of the parameters. So far we presented thermodynamic 
quantities as a function of the temperature parameter r. In Fig. [9] we show 
again Polyakov loop susceptibility xp an d heat capacity C, but now as a 
function of chemical potential fi at fixed r = 0.12 and n = 0.02. In other 
words, in Figs. [3] and 0] this corresponds to a horizontal cut through the 
k = 0.02 crossover region at r = 0.12. From the fact that the maxima of xp 
and C do not coincide in that parameter region we have already concluded 
in the discussion of Fig. H] that this is a crossover region. The absence of 
volume scaling for xp an d C confirms this conclusion. 
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7 Summary 



In this paper we have presented a Monte Carlo simulation of the SU(3) 
spin model in a flux representation. Using the flux degrees of freedom the 
complex phase problem is removed and a Monte Carlo simulation becomes 
possible also at finite chemical potential. After a suitable reparametrization 
of the flux representation |2j we presented a Monte Carlo algorithm that 
alternates updates for the two kinds of monomer and dimer variables. 

Our approach was carefully tested: For /i = 0we compared our results 
to a Monte Carlo simulation using the conventional spin representation. For 
small r and arbitrary \i our Monte Carlo results perfectly match the outcome 
of a perturbative expansion in r. These tests confirm the validity of the map 
to the flux variables, the correct representation of the observables and test 
the implementation of the algorithm. 

We also compared the recently published results from a complex Langevin 
simulation of the SU(3) spin model [7] to the data from our flux simulation. 
We find very good agreement between the two methods which is a valuable 
test for both, the flux and the complex Langevin approach. 

Having tested and compared the implementation of the Monte Carlo 
method in the flux approach we continued with a systematic analysis of the 
phase diagram in the T-fi plane. The phase boundaries were determined from 
the maxima of the Polyakov loop susceptibility xp an d the heat capacity 
C. We found that for small k and [i the transition is of first order. In these 
cases xp an d C display volume scaling, the histograms show a two-state 
signal and the positions of the maxima of xp an d C coincide. For larger 
values of k and /i we only found crossover behavior. The crossover region 
widens with increasing fi, i.e., the positions of the maxima of xp an d C 
become more separated. 

The analysis of the SU(3) spin model in this paper is the first complete 
mapping of the phase diagram for a non-abelian gauge group. We hope that 
the techniques developed here provide a useful step towards new approaches 
for the simulation of QCD or QCD-like systems with chemical potential. 
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